###
### Function to generate alternative roll call matrix
###

   # matching on
   #     + pattern of missingness
   #     + pairwise item correlations
   #     + item response frequencies
   # but using latent multivariate normal response model

max_det_cor_impute <- function(Sigma){
    
    n_vars <- ncol(Sigma)
    
    # identify the number of missing elements of Sigma
    n_missing <- sum(is.na(Sigma))/2
    
    # calculate negative determinant, based on filled in imputes
    detf <- function(imputes){
        Sigma_Imputed <- Sigma
        i <- 1
        for (j in 1:(n_vars-1)){
            for (jj in (j+1):n_vars){
                if (is.na(Sigma[j,jj])) {
                Sigma_Imputed[j,jj] <- Sigma_Imputed[jj,j] <- imputes[i]
                i <- i + 1
                }
            }
        }
        return(-det(Sigma_Imputed))
    }
    
    # use optim to find maximum determinant solution
    start_tmp <- mean(as.vector(Sigma)[as.vector(Sigma) != 1],na.rm=TRUE)
    start_values <- rep(start_tmp,n_missing)
    optim_out <- optim(start_values,detf,method="L-BFGS-B",
                       lower = -1, upper = 1)
    
    # populate Sigma_MD with Max Det solution
    
    Sigma_MD <- Sigma
    i <- 1
    for (j in 1:(n_vars-1)){
            for (jj in (j+1):n_vars){
                if (is.na(Sigma[j,jj])) {
                Sigma_MD[j,jj] <- Sigma_MD[jj,j] <- optim_out$par[i]
                i <- i + 1
                }
            }
    }
    
    # if not positive definite, truncate negative eigenvalues and rescale to valid correlation matrix
    
    n <- dim(var(Sigma_MD))[1L]
    E <- eigen(Sigma_MD)
    Sigma_CM1 <- E$vectors %*% tcrossprod(diag(pmax(E$values, 0), n), E$vectors)
    Balance <- diag(1/sqrt(diag(Sigma_CM1)))
    Sigma_CM2 <- Balance %*% Sigma_CM1 %*% Balance  
    
    # return valid correlation matrix
    
    return(Sigma_CM2)
    
}

polychor_pairwise <- function(mat){
    
    # calculate observed polychoric correlations for items using observed data
    n_vars <- ncol(mat)
    Sigma <- matrix(NA,n_vars,n_vars)
    # populate diagonal
    for (j in 1:n_vars) Sigma[j,j] <- 1
    # calculate pairwise polychoric correlations
    for (j in 1:(n_vars-1)){
        for (jj in (j+1):n_vars){
            if (j != jj) {
                # if pair of items appears, set correlation to observed
                tmp <- suppressWarnings(polychor(mat[,j],mat[,jj]))
                if (!is.na(tmp)) Sigma[j,jj] <- Sigma[jj,j] <- tmp
            }    
        }    
    }
    
    return(Sigma)
}

mvn_rc_generator <- function(mat){
    
    n_vars <- ncol(mat)
    
    # calculate observed polychoric correlations for observed pairs
    Sigma <- polychor_pairwise(mat)
    
    # impute unobserved pairwise correlations to maximise correlation matrix determinant
    # see https://royalsocietypublishing.org/doi/full/10.1098/rsos.172348#FN2R for justification
    Sigma_imputed <- max_det_cor_impute(Sigma)
    
    # generate multivariate normal random latent votes
    y_star <- mvrnorm(nrow(mat),rep(0,n_vars),Sigma_imputed)
    
    # dichotomise latent votes into observed votes
    y <- matrix(NA,nrow(mat),ncol(mat))
    for (j in 1:n_vars) y[,j] <- y_star[,j] > qnorm(prop.table(table(mat[,j]))[1])
    
    # match missingness of original data set
    mat_alt <- y * (mat >= 0)
    
    # return new roll-call matrix
    return(mat_alt)
}

